options(scipen = 999)

# Load the necessary libraries
library(Synth)
library(dplyr)
library(gsynth)
library(writexl)

data <- read.csv("base_merged_PWT_PERCAPITA.csv", row.names = 1) #GDP Per Capita PPP 2017 - PWT 10.01 - Output-side real GDP at chained PPPs (in mil. 2017US$)
data <- data %>% subset(year <= 2009)

data$country_name <- as.character(data$country_name)
data$country_code <- as.character(data$country_code)
# Swapping the 3rd and 4th columns along with their variable names
# Store the variable names
gdp_percapita <- names(data)[3]
country_name <- names(data)[4]
# Swap the variable names
names(data)[3] <- country_name
names(data)[4] <- gdp_percapita
# Swap the columns
temp <- data[, 3]  # Store the 3rd column temporarily
data[, 3] <- data[, 4]  # Replace the 3rd column with the 4th column
data[, 4] <- temp  # Replace the 4th column with the stored 3rd column

country_left = c('China', 'Philippines', 'Canada', 'United States', 'South Africa','Australia') # BENCH

treated_country = 'Chile'
intervention_year = 1990 #1990
# Set the treated unit (country) index for which you want to estimate the treatment effect
# treated_country <- 'Chile'  # Replace this with the index of the country you want to treat

# Identify the pre-treatment and post-treatment periods
pre_treatment_period <- c(1975:(intervention_year-1))   # Replace with the years of your pre-treatment period
post_treatment_period <- c(intervention_year: 2009)  # Replace with the years of your post-treatment period

# Create the outcome and covariates matrices
outcome <- data[data$country_name == treated_country, "gdp_percapita"]

covariate_names <- c(names(data)[5:16])

# Chequeando variables excluidas
names(data)[!names(data) %in% covariate_names]

covariates <- data %>% subset(country_name != treated_country) %>% select(all_of(covariate_names))

data$country_index <- as.factor(data$country_code)
data$country_index <- as.numeric(data$country_index)

# Guarda los resultados

library(car)
model <- lm(gdp_percapita ~ ., data=data %>% select(c(all_of(c('gdp_percapita', covariate_names)))))
vif_values <- vif(model)
print(vif_values)

treated_country_index = unique(data$country_index[which(data$country_name == treated_country)]) #dejar afuera tratado
left_country_index = unique(data$country_index[which(data$country_name == country_left)]) #dejar afuera leave one out
if(length(treated_country_index) > 1){treated_country_index = treated_country_index[1]}

input_dataprep_para_SC <-				
  dataprep(data,				
           predictors = covariate_names,				
           dependent     = "gdp_percapita",				
           unit.variable = "country_index",				
           time.variable = "year",				
           unit.names.variable = "country_name",				
           treatment.identifier  = treated_country_index,				
           controls.identifier   = c(1:24)[!c(1:24) %in% c(treated_country_index, left_country_index)], #all the 98 regions from the donor pool				
           time.predictors.prior = pre_treatment_period,				
           time.optimize.ssr     = pre_treatment_period,				
           time.plot             = c(pre_treatment_period,post_treatment_period))				


modelo.SC <- synth(input_dataprep_para_SC)	

tablas_resultados.SC   <- synth.tab(				
  dataprep.res = input_dataprep_para_SC,				
  synth.res    = modelo.SC)				

# Pesos países
w.at11 <- as.data.frame(tablas_resultados.SC$tab.w)				
w.at11				

w_suficiente <- w.at11 %>% subset(w.weights > 0.001)

# Importancia relativa variables
v.at11 <- as.data.frame(tablas_resultados.SC$tab.v)				
v.at11				

v_suficiente <- v.at11 %>% subset(v.weights > 0.001)


# TEST PLOT! -------------------------------------------------------  

#-----Matching treated vs Synth--------
matching <- as.data.frame(tablas_resultados.SC$tab.pred)	
matching

synth.out.at11 <- synth(input_dataprep_para_SC)	

synth.tables.at11   <- synth.tab(				
  dataprep.res = input_dataprep_para_SC,				
  synth.res    = synth.out.at11)	

# Crear un gráfico para visualizar los resultados --------
path.plot(
  dataprep.res = input_dataprep_para_SC, # los datos preparados
  synth.res = synth.out.at11,       # los resultados de la función synth()
  Main = "Chile's GDP Per Capita",    # el título principal del gráfico
  Ylab = "GDP per capita",         # la etiqueta del eje y
  Xlab = "Year",                    # la etiqueta del eje x
  Legend = c("Chile", "Synthetic Chile"), # la leyenda
  # CI = TRUE                        # si quieres incluir un intervalo de confianza
)


# Agrega una línea vertical en el año 1990 (color rojo)
abline(v = 1990, col = "red") #1990
  
  
